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Abstract 

Monte  Carlo  simulations  are  reported  for  the  solvent  primitive  model  (SPM)  of  the 
electrical  double  layer  at  electrically  charged  and  neutral  surfaces.  Both  solvent  and  ions 
are  modeled  as  hard  spheres  with  interactions  governed  by  hard  sphere  and  Coulomb 
interactions.  Density  profiles  of  solvent  and  ions  and  electrostatic  potential  profiles  are 
presented  for  1:1  electrolyte  concentrations  of  1  and  2M  at  300K.  Compaiision  of  results 
at  charged  and  neutral  walls  indicates  that  the  solvent  structure  contributes  significantly 
to  the  layering  of  ions  and  decrease  of  potentials  in  the  double  layer,  a  result  not  obtainable 
from  simulations  of  the  conventional  primitive  model  (PM)  of  the  double  layer.  The  solvent 
induced  steric  ordering  of  electrolyte  ions  at  a  neutral  wall  can  create  a  space  charge  layer. 
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1.  Introduction 


Separation  of  charge  in  the  vicinity  of  a  surface  that  is  either  charged  or  neutrsd  oc¬ 
curs  in  a  variety  of  phenomena  in  electrochemical,  biological,  tribological,  and  colloidal 
sciences.  For  low  surface  charge  and  dilute  electrolyte  solutions,  the  classical  theory  of 
Gouy  [1]  and  Chapman  [2],  and  later  modified  by  Stem  (GCS)  [3]  gives  predictions  of  the 
interfacial  ion  and  potential  distributions  that  are  in  excellent  agreement  with  experiment 
and  computer  simulations.  More  recent  theoretical  and  computational  studies  have  demon¬ 
strated  that  the  detailed  structure  of  the  solvent  and  ion  are  fimdamental  to  accurately 
describe  the  interfacial  double  layer  properties  for  highly  coupled  systems  at  high  bulk  ion 
densities.  Among  the  models  of  the  electrical  double  layer,  the  most  widely  used  is  the 
primitive  model  (PM),  in  which  ions  are  represented  as  point  charges  that  are  imbeded  in 
hard  spheres,  the  solvent  is  modeled  as  an  isotropic  dielectric  continuum,  and  the  surface 
modeled  as  a  hard  wall  with  uniform  surface  charge  density.  Vigorous  investigation  of  the 
PM  by  computer  simulations  [4*11],  integral  equaticm  methods,  and  functional  theories 
[12-19]  diiring  the  last  decade  have  revealed  phenomena  of  the  electrical  double  layer  that 
are  not  capttired  by  either  the  GC  or  GCS  theories.  For  high  surface  charge  densities  or 
high  electrolyte  concentrations,  the  ion  density  profiles  in  the  PM  predict  a  mudi  more 
highly  organized  layering  of  ions  at  the  charged  surface  than  either  the  GC  or  GCS.  In 
particular,  for  1:1  electrolytes,  the  conterion  profile  shows  two  distinct  layers.  For  both 
1:1  and  2:2  electrolytes  at  high  concentrations,  there  is  charge  inversion  with  coion  den¬ 
sities  exceeding  the  counterion  density  in  the  second  layer.  Neither  of  these  phenomena 
are  predicted  by  the  GCS  model,  although  the  probable  occurrence  of  both  were  discussed 
by  Grahame  [20]  in  the  interpretation  of  the  interfacial  tensicm  measurements  of  Hg.  The 
obvious  disadvantage  of  the  PM  is  that  the  interfadal  structure  of  the  solvent  is  minimized, 
although  it  is  well  known  that  ion  solvation  and  dipole  orientation  at  the  surface,  both  of 
which  are  critical  in  determining  the  potential  and  ion  distributions,  are  correlated  with 
the  molecular  structure  of  the  solvent.  The  development  of  realistic  models  of  water  and 
the  eiifUised  primitive  model  [21],  in  which  the  solvent  is  treated  as  a  dipole  imbedded  in 
a  hard  sphere,  allows  the  effects  of  the  solvent  structure  on  the  double  layer  structure  to 
be  studied.  The  complexity  ct  realistic  models  of  water  limits  studies  of  these  systems  to 
computational  simulations;  however,  very  recent  investigations  of  an  electrified  PtfHiO 
interface  [22]  (no  ions  present)  demonstrates  the  importance  molecular  solvent  structure 
on  the  interface  potential  distribution. 
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A  recent  theoretical  study  of  the  civilized  primitive  model  of  the  electrical  double  layer 
by  Russier  et  al.  [23]  notes  the  importance  of  steric  effects  imposed  by  the  solvent  on  the 
ion  densities  expected  at  a  hard  wall.  These  authors  showed  that  the  closest  approach  of  an 
ion  to  a  neutral  wall  is  determined  by  electrostatic  ion  dipole  interactions  between  ion  and 
solvent  which  is  proportional  to  the  inverse  square  of  the  distance  between  the  interacting 
hard  spheres  and  the  solvent  diameter.  The  results  of  Russier  et  2d.  are  qusmtitatively 
different  from  any  previous  restilts  foimd  in  the  PM,  in  that  the  closest  approach  of  the 
ion  does  not  neccesarily  correspond  to  intimate  contact  between  the  wall  and  ion. 

In  this  work,  we  present  Monte  Carlo  simulations  of  a  solvent  primitive  model  (SPM) 
of  the  electrical  double  layer  in  which  the  effects  due  to  the  finite  size  of  solvent  molecules 
are  included.  Ion  density  and  potential  profiles  are  compared  for  electrolytes  in  contact 
with  charged  and  neutral  walls  for  several  surface  charge  densities  and  electrolyte  concen¬ 
trations.  The  key  findings  of  our  present  work  is  that-  the  finite  size  of  the  solvent  results 
in  highly  ordered  layering  of  ions  at  both  charged  and  neutral  walls.  For  electrolytes  com¬ 
prised  of  different  size  anions  and  cations,  the  solvent  induced  structure  results  in  a  space 
charged  layer  about  two  molecules  thick  at  a  neutral  wall. 

11.  Model  System 


The  systems  simulated  in  this  work  are  mixtures  of  1:1  primitive  electrolytes  and  hard 
sphere  solvent  molecules  contained  between  two  infinite  planar  surfaces.  One  of  the  surfaces 
(at  r  =  0)  is  charged  and  the  other  (at  z  =  H)  is  neutral.  The  systems  are  represented 
computationally  as  a  periodic  array  of  closed  rectangular  boxes  of  height  H  and  a  square 
base  of  length  L  on  a  side.  The  surface  charges  on  the  wall  at  r  s  0  is  balanced  by  the 
excess  counterions  in  the  solution.  N+  and  and  and  are  the  number  and  valence 
of  cations  and  anions  in  the  enclosed  system.  Electroneutrality  requires 

\N.q.  -  (1) 


where  <j  is  the  sur&ce  charge  denaty.  The  interaction  potential  between  particles  i  and  j 


rr/  ^  + 


and  between  particle  i  and  the  wall  is, 


■■)  =  {_ 


00 ,  Zi<  di/2 

2‘K<Tqiezi/€q  Zi>  di/2, 


3 


$ 


0 


where  Zi  is  the  distance  between  the  wall  and  particle  i.  In  equations  (2a,b),  d,  is  the 
diameter  of  particle  i  ,  e  is  the  electronic  charge,  e  is  the  dielectric  constant  (taken  to  be 
78.5  ,  the  value  of  bulk  liquid  water,  in  the  simulation).  Of  course,  if  either  i  or  j  is  a 
neutral  particle  the  Coulomb  potential  in  (2a)  or  (2b)  is  zero. 

In  the  simulations  the  system  is  periodically  extended  in  the  x  and  y  directions.  The 
resulting  long  range  Coulomb  potential  between  ions  and  their  periodic  replicas  is  given 
by  the  sum  obtained  by  Lekner  [24,7] 

oo  oo 

U{rij)  =  4  ^  cos(/^)  ^2  +  2xm)*  +  —  ln(coshr7  -  cosC).  (3) 

/a=l  -00 

with  the  notation  that  (  —  2ir(xi  —  Xj)/L,  rf  =  2x(yi  —  yi)/L,  and  C  =  2n{zi  —  Zj)/L. 
Kq{x)  is  the  modified  Bessel  function  which  decreases  to  zero  very  quickly  as  z  increases. 

m.  Simulation 


In  order  to  run  the  simulation  corresponding  to  fixed  bulk  concentrations,  we  first 
run  the  grand  canonical  Monte  Carlo  (GCMC)  simulations  for  the  primitive  electrolyte 
systems  in  the  following  way  [4-5j.  The  addition  and  deletion  are  attempted  with  equal 
probability  of  about  10  percent  and  the  moves  are  attempted  at  80  percent  probability. 
At  each  addition  or  deletion  a  neutral  pair  of  cations  and  anions  are  added  or  deleted.  If 
the  probability  of  accepting  a  trial  step  from  state  i  to  state  j  is  then 

where  and  Nf^  are  the  number  of  cations  and  anions  in  states  i  and  j  and 

B  »  21n7^  +  ln(n+V'+)  +  ln(n_VL)  (5) 

with  y±  the  mean  ionic  activity  and  rij  is  the  number  density  of  particle  i  with  n,-  s  Ni/Vi. 
Vi  is  the  volume  accessible  to  particle  i.  The  Markov  chain  is  set  up  according  to 

fij  =  mm(l,  fij/fji)  addition. 


fji  =  min(l,  fji/fij)  deletion, 
fij  as  mm(l,  exp{—PUij)  move, 
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where  ^  =  l/feT  with  k  is  the  Boltzman  constant  and  T  is  the  temperature.  In  the 
simulations,  T  is  fixed  at  300  Kelvin.  The  ionic  activity  coefficients  ln7±  =  —0.127  and 
0.271  [5]  are  used  in  the  GCMC  simulations  respectively  for  systems  (a  —  c)  listed  in  Table 
I.  For  systems  d  and  c,  there  are  no  previously  available  data  for  ln7±.  Using  equation  (4), 
for  a  1  molar  1:1  electrolyte,  we  obtained  the  ionic  au:tivities  from  bulk  electrolyte  GCMC 
simulations  to  be  In  7±  =  -0.292  and  -0.44  for  systems  d  and  c,  respectively. 

For  the  confined  systems,  the  average  number  of  cations  and  anions  corresponding 
to  the  right  bulk  concentrations  of  the  systems  are  obtained  from  GCMC  simulations. 
Subsequently,  simulations  of  the  SPM  model  of  the  double  layer  were  carried  out  in  the 
canonical  Monte  Carlo  (CMC)  [3]  simulations  tising  systems  (a  —  e).  The  total  number  of 
solvent  molectdes  and  ions  was  fixed  at 

(JVr+  +  iV_  +  No)d^lHL‘^  -  0.7,  (7) 

where  and  iV_  are  the  average  number  of  free  ions  in  the  confined  PM  fluid.  This  total 
number  density  corresponds  roughly  to  an  aqueous  solution  at  300K.  At  such  high  density 
a  grand  canonical  Monte  Carlo  simulation  is  very  hard  to  carry  out  due  to  the  difficulties 
in  inserting  particles.  The  simulation  process  can  be  divided  into  four  steps  : 

(1)  GCMC  for  bulk  ln7i.  and  bulk  concentration  C 

(2)  GCMC  for  confined  system  with  ln7±  and  C  =►  Average  JV+  +  AT- 

(3)  Add  solvent  =►  SPM,  rui®  '''  0.7 

(4)  CMC  for  SPM 

When  the  electrolyte  concentration  is  1  molar,  the  box  has  a  height  of  H  =  12d  and 
base  of  L  =  4.4d;  reduced  surfisce  charge  densities  of  o*  =  ud^/e  —  0.42  and  0.7,  with 
d  =  4.25  A,  were  used  in  the  simulations.  For  an  e]ectrol}rte  concentration  of  2  molar, 
H  =  14.14d  and  L  —  4.5d,  and  a*  a  0.40.  These  parameters  as  well  as  the  solvent  and 
ion  diameters  are  listed  in  Table  I.  The  long  range  Coulomb  interaction  is  handled  by  the 
summation  technique  given  in  section  11.  The  probability  of  acceptance  used  in  the  SPM 
simulations  is  10  percent.  During  the  run,  the  systems  were  equilibrated  for  50,000  steps 
and  then  was  run  (0.5  —  1)  x  10*  to  accumulate  averages  to  calculate  the  density  profiles 
n(z)  and  the  mean  electroetatk  potential  ^(z).  The  density  n(z)  was  obtained  counting 
the  number  of  particles  in  slices  parallel  to  the  charged  wall  and  dividing  by  the  volume. 
Very  long  runs  are  required  compared  to  the  primitive  electrolytes  simulations,  in  <^der  to 
eliminate  the  noise  in  the  density  profiles.  0(z)  is  obtained  eta.  the  relation 
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Here  we  take  the  right  wall  as  the  reference  point  of  the  potential. 


IV.  Results  and  Discussion 


In  the  following  we  report  our  simvilation  results  of  the  density  profiles  emd  mean 
electrostatic  potentials  for  the  different  systems  that  are  listed  in  table  I.  For  all  the 
systems  studied,  the  left  side  wall  is  charged  and  the  right  wall  (at  z  =  ff  )  is  neutrad.  The 
surface  charge  density,  the  ion  and  solvent  density,  and  the  mean  electrostatic  potential 
are  reported  in  the  dimensionless  forms 

trtfi  n- 

- - ,  n*  =  — ,  and  (9) 

e  n,o 

where  Uio  is  the  bulk  number  density  for  species  i  ion  and  is  the  number  density  for  tne 
neutral  solvent. 

To  compare  the  results  obtained  here  with  those  obtained  for  the  primitive  electrolytes, 
in  figure  1,  we  show  the  SPM  results  as  continuous,  dotted  or  dashed  curves  and  the  PM 
results  as  open  and  filled  drdes  for  the  case  of  o*  s  0.42  and  C  ^  IM.  These  PM  resiilts 
can  be  obtained  in  step  (2)  of  the  simulations  or  from  previous  simulations  [4,7].  As  is 
well  known,  the  Gouy-Chapman  theory  and  simiilations  of  the  primitive  electrolytes  are 
in  i^eement  for  low  surface  charge  density  and  concentration.  Both  the  density  profile 
and  the  mean  electrostatic  potential  are  monotonic.  FVom  the  figure  we  can  see  that,  the 
solvent  structure  included  in  the  SPM  simulation  has  a  dramatic  effect  on  the  electrolyte 
ion  profiles.  The  presence  of  solvent  molecules  at  the  walls  induces  strong  structure  in 
the  i<m  distributions  at  both  charged  and  neutral  walls.  Near  the  left  wall,  the  neutral 
particles  and  counterions  have  five  layers,  the  coions  show  four  layers  and  near  the  right 
wall  there  are  four  layers  for  all  particles.  A  small  residue  of  the  coions  near  the  left  wall  is 
in  the  SPM  model  simulation.  These  structural  features  are  absent  in  the  PM  simulations. 
Figure  2  shows  the  mean  electrostatic  potential  for  the  same  system  as  in  figure  1,  the 
potential  decreases  rapidly  nrar  the  charged  wall  and  decays  to  zero  firom  a  negative  value. 
As  mentioned  in  section  III,  the  noise  in  the  density  profiles  prevents  us  from  getting  more 
accurate  results  to  see  a  clear  trend  of  the  potential. 

Figures  3  and  4  are  firom  the  simulations  with  tr*  =  0.7  and  C  =  lAf .  The  layering 
of  ions  is  very  distinct  as  seen  in  figure  1.  Because  the  strong  field  on  the  left  wall,  the 
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coions  were  completely  repelled  from  the  wall  in  the  PM.  However,  the  hard  sphere  solvent 
and  ion  interactions  in  the  SPM  again  restilt  in  a  significant  ordering  of  coions  near  the 
chEirged  sxuface.  The  mean  electrostatic  potential  of  the  SPM  decreases  more  rapidly  near 
the  charged  wall  than  does  that  of  the  PM. 

In  figures  5  and'fi,  we  change  the  ion  concentration  to  2M  and  use  a*  ~  0.4.  Again  the 
layering  of  ions  is  very  strong  and  similar  to  figures  1  and  2.  Furthermore,  a  charge  inversion 
occured  in  the  second  and  third  layer  with  the  coion  density  exceeding  the  counterion. 

To  show  that  the  highly  structured  ion  distributions  of  the  SPM  are  not  due  to  the  fact 
that  the  solvent  and  charged  particles  have  the  same  diameter,  we  ran  simulations  for  two 
different  electrolytes  comprised  of  dissimilar  sized  particles  d+  =  0.75d_  and  d+  =  0.5d_  . 
The  results  are  given  in  figures  7  through  10.  These  results  show  that  the  coion  peaks  shift 
due  to  the  size  effect  but  there  is  still  significant  layering  structure  in  the  ion  and  solvent 
density  distributions.  Ion-size  induced  charge  separations  is  clearly  seen  near  the  neutral 
wall.  Due  to  the  noise  in  the  density  profiles,  we  cannot  make  an  accurate  estimate  of 
the  potentials  resulting  from  these  separations.  Comparing  the  density  profiles  in  figures 
7  and  9,  one  can  see  that  the  layering  seems  miunly  determined  by  the  larger  solvent  and 
counterions  near  the  two  walls,  even  though  the  effect  of  the  coion  size  made  the  peak  shift 
a  little.  Again  these  structures  are  absent  in  PM  electrolytes. 

FVom  the  above  analysis  and  comparision  with  the  primitive  electrolyte  simulations 
we  see  the  essential  role  of  the  solvent  in  determining  the  structure  of  electrolytes  at  solid 
surfaces.  In  fact,  the  layering  structures  appearing  in  the  SPM  electrolytes  are  mainly  due 
to  the  presence  of  the  solvent  molecules.  The  small  residue  of  coions  near  the  charged 
surfaces  is  due  to  the  layering  of  the  solvent  and  counterions.  The  fact  that  we  did  not 
observe  stronger  structvires  near  the  charged  wall  at  higher  surface  charge  density  (<r*  = 
0.7)  than  at  lower  sxirface  charge  densities  {  a*  0.4)  confirms  the  idea  that  the  effect  of 
the  solvent  plays  a  more  important  role  in  determining  the  electrolyte  structure  than  solely 
by  the  surface  charge  hard  sphere  ion  interactions,  as  is  the  case  of  the  primitive  electrolyte. 
In  all  cases  that  have  been  studied  in  this  paper,  the  mean  electrostatic  potential  of  the 
SPM  is  smaller  than  that  obtained  from  the  PM  near  the  charged  wall. 

The  apparent  structure  of  the  interfaces  obtained  from  the  solvent  primitive  model 
may  provide  new  insight  into  the  interaction  of  charged  particles  near  surfaces.  The  layering 
of  molecules  near  a  sxirface  would  be  reflected  by  oscillatory  interaction  forces.  Force 
measurements  using  mica  surfaces  have  revealed  short  range  oscillatory  electrical  double 
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layer  forces  [25-29]. 
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Table  I.  Parameters  used  in  this  paper,  a*  =  a<P ft,  side  length  i,  surface  separation 
H.  Number  of  cations,  anions,  and  neutral  particles  are  Nj^  ,  iV_,  and  No,  and  Nc  = 
N+  +  N_,  AN  =  N_  —  N+  =  aA^  A  =  L  x  L.  C  is  the  bulk  ion  concentration  for  PM 
electrolytes  in  molar.  The  number  of  steps  run  in  the  simulations  are  given  as  Nsteps  in 
unit  of  10®  steps  in  the  last  column.  In  all  the  cases  studied  in  -this  paper  the  neutral 
solvent  particle  and  the  anion  have  diameter  of  do  =  «i_  =  d  =  4.25  A.  For  systems  a-b,  c, 
d  and  e,  ln7±  =  —0.127,  0.271,-0.292  and  —0.44,  respectively. 


sys 

H/d 

L/d 

d+/d. 

No 

No 

AN 

C{M) 

Nsteps/ 10 

a 

4.36 

26 

134 

8 

1.0 

b 

4.47 

H 

32 

136 

14 

0.5 

c 

14.14 

4.5 

0.40 

52 

148 

8 

2.0 

0.6 

d 

12.0 

4.36 

0.42 

26 

134 

8 

1.0 

0.5 

e 

12.0 

4.36 

0.42 

26 

134 

8 

1.0 

0.5 

Figure  Captions 


Figure  1. 


Figure  2. 


Figure  3. 


Figure  4. 


Figure  5. 


Figure  6. 


Figure  7. 


Figure  8. 


Figure  9. 


Figure  10. 


Reduced  Density  profiles  n*{2)  =  n{z)lnQ  for  a  1:1  electrolyte  at  C=1M  and  <i+  = 
d-  =  d  =  4.25  A.  The  left  wall  has  a  surface  charge  density  of  a*  =  ad?/e  =  0.42 
zind  the  right  wall  is  uncharged.  Solvent  primitive  model  (SPM)  simulations:  solvent 

( - )  counterions  ( — ),  coions  (...).  Primitive  model  (PM)  simulations:  open  (o) 

and  filled  (•)  circles  correspond  to  counter  and  coions,  respectively. 

Reduce  Mean  electrostatic  potential  0*(z)  for  a  1:1  electrolyte  for  the  same  conditions 
as  in  figure  1.  SPM  simulations  ( - )  ;  PM  simulations  ( — ). 

Reduced  Density  profiles  n*{z)  for  a  1:1  electrolyte  at  C=1M  and  a*  =  0.7  .  Other 
conditions  and  keys  are  the  same  as  in  figure  1. 

Reduce  Mean  electrostatic  potential  for  a  1:1  electrolyte  for  the  same  conditions 
as  in  figure  3.  SPM  simulations  ( - )  ;  PM  simulations  ( — ). 

Reduced  Density  profiles  n*{z)  for  a  1:1  electrolyte  at  C=2M  and  a*  =  0.40  .  Other 
conditions  and  keys  are  the  same  as  in  figme  1. 

Reduce  Mean  electrostatic  potential  0*(z)  for  a  1:1  electrolyte  for  the  same  conditions 
as  in  figure  5.  SPM  simulations  ( - )  ;  PM  simiilations  ( — ). 

Reduced  Density  profiles  n*(r)  for  a  1:1  electrolyte  at  C=1M  and  a*  =  0.42  .  d^.  = 
|d_.  Other  conditions  and  keys  are  the  same  as  in  figure  1. 

Reduce  Mean  electrostatic  potential  ^*{z)  for  a  1:1  electrolyte  for  the  same  conditions 
as  in  figure  9.  SPM  simulations  ( - )  ;  PM  simulations  ( — ). 

Reduced  Density  profiles  n*{z)  for  a  1:1  electrolyte  at  C=1M  and  a*  =  0.42  .  d+  = 
^d_.  Other  conditions  and  keys  are  the  same  as  in  figure  1. 

Reduce  Mean  electrostatic  potential  ^*(r)  for  a  1:1  electrolyte  for  the  same  conditions 
as  in  figure  7.  SPM  simtilations  ( - )  ;  PM  simulations  ( — ). 
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Mean  Electrostatic  Potential  Density  Distribution  n  (z) 
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Ilean  Electrostatic  Potential  Density  Distribution  n*(z) 


Mean  Electrostatic  Potential  Density  Distribution  n*(z) 


